Arbitrarily High Order BGK-Shakhov Method for the Simulation of Micro-Channel Flows

A new arbitrarily high order method for the solution of the model Boltzmann equation for micro-channel flows in the transitional regime is presented. The Bhattnagar-Gross-Krook approximation of the Boltzmann collision integral is implemented, with Shakhov’s modification, and the resulting system of equations solved by a discrete ordinate method. The method approximates velocity space using a truncated Hermite polynomial expansion of arbitrary order and performs the associated integration by Gauss-Hermite quadrature. This approach conserves mass, momentum and energy during relaxation of the discretised velocity space towards equilibrium. Physical space is discretised by discontinuous Legendre polynomial expansions with both the spatial representation and conservative flux calculation being of arbitrary order. Owing to the high order spatial representation of the discretised velocity space the BGKShakhov relaxation process is carried out in a ‘continuous in space’ manner. New high order boundary conditions of the inviscid slip wall and no-slip wall are implemented. A new fully diffuse reflection boundary condition, built on the high order spatial information available in the method, is also proposed. Results are presented for low speed planar Couette flow and non-linear channel flow. Introduction The Boltzmann equation, which describes the evolution of the velocity distribution of a dilute monatomic gas through binary elastic collisions, allows the simulation of flows through the transitional regime, and beyond, where the Navier-Stokes equations are not valid. Methods that allow the efficient solution of the Boltzmann equation are therefore highly sought after. The Boltzmann equation, ∂ f ∂t +~ξ · ∂ f ∂~x = Ω( f )≈ ν ( f target− f ) , (1) describes the gas flow in terms of the velocity distribution function, f = f ( ~x,~ξ, t ) , which is a function of position, ~x, microscopic absolute velocity, ~ξ, and time, t. The non-linear integral collision term, Ω( f ), is numerically expensive to calculate and can be replaced with the much simpler Bhattnagar-GrossKrook (BGK) model [1], as shown, which relaxes f towards a known target distribution, f target, with relaxation frequency ν. To circumvent the fixed Prandtl number (Pr = 1) limitation of the BGK model, the Shakhov target distribution may be used [2]. The Shakhov distribution includes the heat flux vector, ~q, as an input which is itself a moment of f . To solve Eq. (1) numerically, the discrete ordinate method (DOM) [3] may be used, whereby the continuous velocity distribution function is replaced by multiple discrete distributions, fi (~x, t), each with a corresponding constant advection velocity, ~ξi. To enforce conservation of mass, momentum and energy, in an efficient manner, the continuous distribution function, f ( ~x,~ξ, t ) , may be approximated as a truncated Hermite polynomial [4]. The discrete distributions of the DOM then correspond to the prescribed abscissa of the Gauss-Hermite quadrature rule selected to allow exact integration of the Hermite polynomial. The linear advection component of Eq. (1) may be implemented using any number of techniques including finite difference, volume, and element approaches. Recently, a scheme that uses discontinuous truncated Legendre polynomials to represent high order spatial variations has been proposed [5]. This method uses the high order information within each finite volume cell to perform linear advection and so presents a uniformly high order scheme that is well suited to the DOM approach. In this paper we give a brief description of the numerical method before outlining the high order boundary conditions that have been implemented and their related issues. Results for channel flows are then displayed with an investigation into the order of convergence of the method and the effect of rarefaction on heat flux into the wall for planar Couette flow. The Numerical Method In this section a brief overview of the numerical method will be presented, for a more complete description refer to Bond et al. [6]. The Boltzmann equation, given in Eq. (1), can be written in discrete form while maintaining exact recovery of all moments of the velocity distribution function up to some arbitrary order, by mapping f ( ~x,~ξ, t ) onto a Hermite subspace according to the method of Shan et al.[4]. The continuous form of the equations is then reduced to a set of discrete distributions, one for each of the velocity abscissa,ξk, required by the Gauss-Hermite quadrature rule. By replacing the collision operator, Ω, with the BGK-Shakhov approximation and a reduction of the spatial dimension, by the method of Chu [7], the Boltzmann equation can be written according to Eqs. (2). The g and h terms relate to the translational and thermal energy components, respectively, of the original distribution function, f , while the gk and h S k denote the use of the Shakhov model as the target distribution. ∂gk ∂t +~ξ ∂gk ∂~x =ν ( gk −gk ) ∂hk ∂t +~ξ ∂hk ∂~x =ν ( hk −hk ) (2) From this point on it should be noted that all vector quantities refer to only two spatial variables i.e.~ξ = [u,v]. Nondimensionalisation We define a characteristic length, L, and speed, C∞, and the characteristic time, t∞, follows accordingly in Eq. (3). A characteristic density, ρ∞, and temperature, T∞ are also required. C∞ = √ RT∞ t∞ = L C∞ (3) The non-dimensional variables for absolute molecular velocity, ~ξ, mean macroscopic velocity, ~Ξ, relaxation frequency, ν, and the stress,~τ, heat flux, ~q, and acceleration, ~a, vectors can then be introduced in Eqs. (4). ~̂x =~x/L ~̂ξ =~ξ/C∞ ~̂Ξ =~Ξ/C∞ t̂ = t/t∞ ν̂ = νt∞ ρ̂ = ρ/ρ∞ θ̂ = T/T∞ ~̂τ =~τ/ρ∞C ∞ ~̂q =~q/ρ∞C 3 ∞ ~̂a =~a/(L/RT∞) ĝ = g/ ( ρ∞/C ∞ ) ĥ = h/ρ∞ Kn = λ/L (4) The Knudsen number, Kn, indicates the degree of translational non-equilibrium in the flow. All quantities from this point on, unless otherwise specified, will be in non-dimensional form and the carats will be omitted. The Shakhov distributions for gk and hk, used as the target distribution in Eqs. (2), can then be given in non-dimensional form according to Eqs. (5-6). gk = g M k [ 1− (Pr−1)(~q ·~ck)(~ck ·~ck−4θ) 5θ3ρ ] (5) hk = g M k θ [ 1− (Pr−1)(~q ·~ck)(~ck ·~ck−2θ) 5θ3ρ ] (6) gk = ρ 2πθ exp ( −k ·~ck 2θ ) (7) The peculiar velocity,~ck, is given by~ck =ξk−~Ξk and the nondimensional collision frequency, Eq. (8), which determines the rate of relaxation can also be defined. The variable ω is the exponent of the power law viscosity although other viscosity laws can be used. ν = ρ Kn √ π 2 θ1−ω (8) Spatial Representation The solution of the set of discrete equations given in Eqs. (2), by the process of linear advection and relaxation, allows for the implementation of an explicit linear advection solver. The Conservative Flux Approximation (CFA) method of Latorre et al.[5] solves the linear advection equation in multi-dimensional situations in a high order manner. The CFA method representats each discrete distribution, gk and hk, as a truncated expansion in Legendre polynomials, of order Nl , in space that is discontinuous across cell interfaces. This can be expressed for the k-th distribution in cell (i, j) according to Eqs. (9-10). gi, j,k (~x) = Nl ∑ m=1 Nl ∑ n=1 a i, j,k Lm (x)Ln (y)≈ g exact i, j,k (~x) (9) a i, j,k = (2m−1)(2n−1) 4 ∫∫ 1 −1 gexact i, j,k (~x)Lm (x)Ln (y)dxdy (10) The polynomial representation given by Eq. (9) is encoded in the coefficient matrix a, of which a(m,n) is the coefficient corresponding to the product of the mth and nth order Legendre polynomials. The exact expression for the advection of the polynomial representaton is given by Eq. (11) and can be seen, for a 1D example with positive advection velocity, in the first two panels of Fig. 1. gk (~x, t +∆t) = gk ( ~x−ξk∆t, t ) (11)